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Quantum vs. Geometric Disorder in a Two-Dimensional Heisenberg Antiferromagnet 
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We present a numerical study of the spin- 1/2 bilayer Heisenberg antiferromagnet with random 
inter-layer dimer dilution. From the temperature dependence of the uniform susceptibility and a 
scaling analysis of the spin correlation length we deduce the ground state phase diagram as a function 
of nonmagnetic impurity concentration p and bilayer coupling g. At the site percolation threshold, 
there exists a multicritical point at small but nonzero bilayer coupling g m = 0.15(3). The magnetic 
properties of the single-layer material La2Cui_ p (Zn,Mg) p 04 near the percolation threshold appear 
to be controlled by the proximity to this new quantum critical point. 

PACS numbers: 75.10. Jm,75.10.Nr,75.40.Cx,75.40.Mg 
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Quantum phase transitions (QPTs) in the presence 
of disorder are the subject of considerable current in- 
terest as they exhibit rich new physics, but are rather 
poorly understood. Experimental examples include the 
cuprate superconductors [Q, heavy fermion compounds 
metal-insulator J3| and superconductor- insulator 
transitions quantum Hall effect and quantum 
magnets ||. A significant amount of theoretical and 
numerical work has been devoted to the study of the 
random Ising chain in a transverse field , the simplest 
quantum many-body system with quenched disorder, and 
to its analog in higher dimensions || |9) . Partly because 
of its relevance to cuprate superconductivity, the spin- 
1/2 square-lattice Heisenberg antiferromagnet (SLHAF) 
has attracted enormous interest |l0j . While the effects of 
a single nonmagnetic impurity are well understood [jnj, 
there exist few theoretical results for finite impurity con- 
centrations Jl^ |. Given the considerable challenges, the- 
ory for quantum systems with disorder is often guided by 
insight provided by numerical work. 

In the absence of disorder, the nearest-neighbor (NN) 
SLHAF can be driven through a QPT by introducing a 
parameter analogous to the transverse field for the Ising 
model. For example, this is achieved by introducing frus- 
trating next-NN couplings |l3| or by coupling two square 
lattices to form an antiferromagnetic bilayer It had 
been argued that random site or bond dilution of the SL- 
HAF in the extreme quantum limit of spin-1/2 may lead 
to a non-trivial QPT j^, |l5|, |l6 ]. However, recent experi- 
mental and numerical ]17|, 18| work suggests that the 
ground state remains ordered for non-magnetic impurity 
concentrations p up to the site percolation threshold , 
and that the critical cluster at p = p* appears to have a 
nonzero staggered moment [fLsl , which would imply that 
the percolation transition is classical. 

In this Letter, we present numerical results for the site- 
diluted spin-1/2 NN bilayer antiferromagnet, with disor- 
der that is fully correlated between the layers ("dimer" 



dilution). At zero bilayer coupling, this problem reduces 
to the previously studied diluted spin-1/2 SLHAF. By 
increasing the strength of the bilayer coupling from zero, 
we are able to increase quantum fluctuations beyond 
those for the spin-1/2 SLHAF. Our finite-temperature 
results allow us to extract the ground state phase dia- 
gram as a function of the strength of quantum fluctua- 
tions and the degree of disorder, and they reveal a new 
multicritical point at p = p* for small, but nonzero bi- 
layer coupling. This phase diagram resembles that of 
the diluted two-dimensional Ising ferromagnet in a trans- 
verse field ||. The present results are relevant to re- 
cent experimental findings for the model spin-1/2 SLHAF 
La 2 Cu 1 _p(Zn,Mg) p 4 fy. 
We study the Hamiltonian 
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where J (Jj_) is the antiferromagnetic planar (bilayer) 
coupling, the first sum is over all planar NN pairs, and 
(-i = (ej = 1) with probability p (1 — p). We define 
the reduced bilayer coupling g = Jj_/ J, and work with 
units in which J = a = fcs = 37-tB = = 1 (a is the 
lattice constant). If a site is removed on one layer, the 
corresponding site in the adjacent layer is also removed. 
This constraint preserves the percolation properties of 
the square lattice. It also ensures that there are no un- 
paired "dangling" spins. Unlike for the diluted SLHAF 
(g = 0), the uniform susceptibility Xu f° r the bilayer 
then contains no Curie-like term. The temperature de- 
pendence of Xu can then be used to determine the crit- 
ical bilayer coupling g c {p) for the diluted system. We 
use the loop-cluster Monte Carlo method ]n| , which has 
provided good results for the diluted spin-1/2 ladder ]2Cf| 
and SLHAF jl7j] . Simulations are performed on large lat- 
tices of up to JV = 512 x 512 x 2 sites, to temperatures 
as low as T = 1/300, with Trotter numbers of 20/T, and 
by averaging 10 to 100 random configurations. By keep- 
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ing the planar lattice size significantly larger than the 
correlation length £(g,p, T), we are able to avoid finite- 
size effects. Several concentrations are chosen, spanning 
the range from the pure system up to p — 0.5, beyond 
the site percolation threshold p„ = 0.40725379(13) |§l). 
In addition to x u (<7,p,T) and £(<7,p,T), we also com- 
pute the staggered susceptibility Xst(g,P,T). The sus- 
ceptibilities are defined as NT Xu = [((E^f) 2 ) ] and 

NT X st = [((EiC- 1 )' s i) 2 )]> where the sum is over a11 
sites N and [...] indicates the disorder average. The in- 
stantaneous staggered correlation function is given by 
NC(r) = sign{v) [(£ 4 S?S? +r >] , where sign{v) = 1 (- 
1) for r separating sites on the same (different) sublat- 
tice. The correlation length is obtained from the behavior 
C(r) ~ e~' r l/£ at large distances. 

At the critical coupling of the pure bilayer, 
Xm(<?c;0, T) ~ T @. In the quantum disordered phase 
(g > g c ), the strong bilayer coupling leads to a nonzero 
singlet-triplet gap, which prevents long-range order in 
the ground state, and X u is exponentially small at low 
temperature. In the ordered phase (g < g c ), Xu ap- 
proaches a nonzero value proportional to the spin stiff- 
ness p2[ . These different characteristic behaviors have 
allowed the determination of the critical bilayer coupling 
for the pure system from finite- T numerics [ pd] ]. This 
should also be possible at nonzero p, as long as there 
exists an ordered phase. Figure [l] shows some of our 
results for Xu f° r the pure system, at the intermediate 
concentration of p = 25%, and at a concentration just 
below the percolation threshold. Note that T is scaled 
by y, so that results for different concentrations fall into 
approximately the same horizontal range. Fitting the 
low-T data to X u = ci + C2T A gives excellent results, as 
shown by the dashed lines in Fig. [|. We interpret pos- 
itive (negative) values of c\ as indicative of an ordered 
(quantum disordered) ground state. For p = 0, we find 
g c k 2.525, in very good agreement with g c = 2.525(2) 
obtained previously |14| . The critical coupling decreases 
with increasing p, and A decreases from its value of 1 
for the pure system to around 0.7 near the percolation 
threshold. We conclude that order persists at g — 0.10 
even at p = 417/1024 (40.722...%), very close to p = p., 
(40.725...%). 

In the ordered phase, the low-T correlation length is 
expected to be exponential in inverse temperature as for 
the pure SLHAF E while £ - J 7 - 1 /* at the quan- 
tum critical point fl23|, where z is the dynamic critical 
exponent. The crossover temperature Tdev at which £ de- 
viates from a power law scales with the deviation from the 
critical point as Tdev ~ \g — 9c\^ (p can be substituted for 
g for cuts at fixed g). Similarly, in the disordered phase, 
the low-T correlation length approaches a constant value, 
but follows power-law behavior at higher temperatures. 
Figure 2 shows scaling plots for four cuts across the phase 
boundary. For the pure bilayer, Fig. 2(a), for which the 
low-T physics can be mapped to the quantum nonlinear 
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FIG. 1: Temperature dependence of the uniform suscepti- 
bility, (a) Pure bilayer: g c » 2.525 and X«(5c,0, T) ~ T, 
as found previously [Q. (b) p = 0.25: g c ~ 1.42 and 
\u ~ T 0,73 ' 5 '. (c) Just below the percolation threshold: <7 C is 
significantly reduced, and y„ ~ yO- 7 o(7)^ 



cr model [ ^3| , the Euclidean time direction is equivalent 
to the spatial dimensions (Lorentz invariance) and z = 1 . 
At zero temperature, Euclidean time extends to infinity, 
so that the critical exponents of the pure system are those 
of the three-dimensional classical Heisenberg model, for 
which v fts 0.705 and 7 w 1.39 {24j. Since v = <fi/z, we 
therefore have fixed <fi = 0.705 in Fig. 2(a). We obtain 
excellent scaling with g c = 2.5215(10), slightly lower than 
the value obtained from Fig. 1(a) and in previous work 
|l4j| . A similar result is obtained for Xst (not shown) 
using 7 = 1.39. 

At p = 25% [Fig. 2(b)], we find z = 1.07(2) and (j) = 
0.95(5), as well as g c (p = 0.25) = 1.412(6), consistent 
with Fig. 1(b). Scaling as a function of p at fixed g = 1.42 
(not shown) gives good results using the same exponents 
and p c (g = 1.42) = 24.85(15)%. Near p = p*, scaling 
as a function of p is easier to obtain, and our results for 
g = 0.20 and 5 = (single layer) are shown in Figs. 2 
(c) and (d), respectively. For g — 0.20, the exponents 
are z = 1.3(1), <f> = 1.45(10), and p c = 39.9(3)%. Recent 
neutron scattering results for La2Cui_ p (Zn,Mg) p 04 
as well as numerical work jl6[ |l7], |l^| indicate that p c = 
p* for the spin-1/2 SLHAF. Therefore, at g = 0, we fixed 
p c — p*. As is evident from Fig. 2(d), this gives excellent 
scaling, with z = 1.65(5) and <f> = 1.8(1). 

Sandvik (fl| concludes that the percolating cluster of 
the NN SLHAF at p = p* has a nonzero ordered moment 
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FIG. 2: Scaling plots for the correlation length (a) for the 
pure spin-1/2 bilayer Heisenberg antiferromagnet, (b) at fixed 
nonmagnetic concentration p = 0.25, and at fixed bilayer cou- 
plings (c) g = 0.20 and (d) g = (SLHAF). 



per spin, so that even in the extreme quantum limit of 
spin-1/2 the percolation transition at T = is still clas- 
sical. If this picture is correct, then one would expect 
the moment per spin on the critical cluster to vanish at 
some nonzero value g m — g c {p*)- This is indeed consis- 
tent with our findings. Figure 3 shows the extrapolated 
ground state phase diagram, and we obtain g m = 0.15(3). 
The exponents for g — 0.10 are very similar to those at 
g = 0.20, and we estimate z m = 1.33(7), v m = 1.11(14), 
and 7 m = 2.69(31). Calculations of the temperature de- 
pendence of the correlation length at this multicritical 
point agree with this estimate for z m . 

The Harris criterion is often used to determine whether 
the exponents at classical phase transitions should change 
in the presence of weak disorder: for disorder to be rel- 
evant, v of the pure system has to satisfy v < 2/d p5|. 
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FIG. 3: Phase diagram of the diluted bilayer. The dashed line 
is a guide to the eye. The percolation threshold is marked by 
the dotted vertical line. The inset is a close-up view of the 
region near the percolation threshold. 



More generally, it is expected that the disordered sys- 
tem satisfies v > 2/d p6[ . For a quantum phase transi- 
tion, the dimension to use is the number of dimensions 
in which there is disorder, in our case d = 2. Since for 
the pure bilayer v « 0.705 < 1, we expect disorder to be 
relevant |23[ . 

The meaning of the exponents obtained along the 
phase boundary (0 < p < p*) is not entirely clear. As 
argued above, disorder is expected to be a relevant per- 
turbation of the pure system. This should lead to a 
unique set of critical exponents for the randomly-diluted 
bilayer Heisenberg antiferromagnet, different from those 
at g c (0) and g mi and from the classical percolation ex- 
ponents. There are no special values of p between p = 
and p — p*, and the exponents appear to vary contin- 
uously along the phase boundary. This is qualitatively 
similar to the two-dimensional Ising magnet, for which 
quenched disorder is a relevant perturbation, and which 
appears to exhibit (finite-temperature) transitions with 
concentration-dependent exponents p7| . However, these 
Ising exponents also show a slight temperature depen- 
dence, and it has been argued that this indicates that 
these are not the true exponents, which should only be 
observable at temperatures asymptotically close to the 
transition temperature, but rather effective exponents 
p7f . A similar masking might occur in the present sys- 
tem, so that true critical exponents would only be re- 
vealed at asymptotically low temperatures. The criti- 
cal points at g c (0) and g m influence the intermediate- T 
physics of systems nearby in parameter space. One might 
expect that at p w 25% the extracted exponents are close 
to the true critical exponents: p — 25% appears to be far 
away from both the pure system and from the percolation 
threshold, so that disorder effects might be large, while 
£ at the temperatures of our study is still much larger 
than the percolation length, which diverges at p = p*. 
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However, we obtain 
criterion v > 1 26 



i> = <fi/ z = 0.89, in violation of the 
Consequently, we believe that the 
exponents extracted for < p < p, are effective, rather 
than true critical exponents. 

The scaling dimension for a quantum phase transition 
is d q = d+ z, where z can differ from 1 in the disordered 
case. Using z m = 1.33, v m = 1.11, and 7 m = 2.69, the 
hyperscaling relationship (3 — {d q v — ^f)/2 gives (3 m « 0.5 
for the multicritical point. Our corresponding effective 
values for g = are z = 1.65(5), v = 1.09(9), and 
7 = 2.65(27), and /3 « 0.67(13). Previous finite-size scal- 
ing results for the NN SLHAF resulted in spin-dependent 
critical exponents, with z = 2.54(8), v = 1.23(16) and 
j3 = 0.50(7) for spin-1/2 Jig] . Based on our results, 
these exponents have to be viewed as effective exponents 
due to the influence of the nearby multicritical point. 
This is consistent with the claim that the spin-1/2 NN 
SLHAF at p = p* should exhibit asymptotic percola- 



tion critical behavior |18|, for which u 



- 4/3 w 1.33, 

7 P = 43/18 w 2.39, and /3 P = 5/36 « 0.14 @. Exper- 
iments for La2Cui_p(Zn,Mg) p 04 near p — p* revealed 
£ t~" t , with vt ~ 0.7, which is very close to the ef- 
fective value 1/z w 0.61 (g — 0) and to l/z m k 0.75 
{g m = 0.15). 

In summary, we have mapped out the phase diagram 
of the spin-1/2 bilayer Heisenberg antiferromagnet with 
quenched disorder in form of interlayer dimer dilution. 
Varying the bilayer coupling has allowed us to further 
increase quantum fluctuations beyond those of the spin- 
1/2 NN SLHAF and to investigate the joint effects of 
quantum fluctuations and quenched disorder. Our re- 
sults for the pure bilayer agree with previous work. We 
find that the critical coupling g c (p) decreases with in- 
creasing p, but remains nonzero even at the percolation 
threshold p = p*. The point g rn = g c (p*) is a new mul- 
ticritical point, and we obtain estimates of several crit- 
ical exponents. Quenched disorder is expected to be a 
relevant perturbation to the pure quantum system and 
to lead to new critical behavior at nonzero p below the 
percolation threshold. The critical exponents along the 
phase boundary (0 < p < p») appear to change contin- 
uously, but it is likely that the true critical behavior is 
masked by finite-temperature effects, and by the effec- 
tive proximity to either the pure fixed point or the new 
multicritical point. The small value of g m « 0.15 indi- 
cates that the intermediate-temperature properties of the 
spin-1/2 SLHAF La 2 Cui_p(Zn,Mg)p0 4 near p = p* 
may be controlled by this nearby quantum critical point. 
We note, that the experimental system might be best 
described by considering a nonzero, frustrating next-NN 
exchange of about 0.05 — 0.10 J jL7], |29|], which, in effect, 
would place it even closer to the multicritical point. The 
phase diagram of Fig. 3 is qualitatively similar to that 
for the randomly diluted 2D Ising model in a transverse 
field U . We note that recent numerical work , which 
uses a different approach, gives g m (j>*) — 0.16(1) with 



z m = 1.28(2), in very good agreement with our results. 
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